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Abstract. The investigation of phase coexistence in systems with muhi-component 
order parameters in finite systems is discussed, and as a generic example, Monte 
Carlo simulations of the two-dimensional q-state Potts model (q=30) on L x L 
square lattices (40 < L < 100) are presented. It is shown that the microcanonical 
ensemble is well-suited both to find the precise location of the first order phase 
transition and to obtain an accurate estimate for the interfacial free energy between 
coexisting ordered and disordered phases. For this purpose, a microcanonical version 
of the heatbath algorithm is implemented. The finite size behaviour of the loop 
in the curve describing the inverse temperature versus energy density is discussed, 
emphasizing that the extrema do not have the meaning of van der Waals-like " spinodal 
points" separating metastable from unstable states, but rather describe the onset of 
heterophase states: droplet/bubble evaporation/condensation transitions. Thus all 
parts of these loops, including the parts that correspond to a negative specific heat, 
describe phase coexistence in full thermal equilibrium. However, the estimates for 
the curvature-dependent interface tension of the droplets and bubbles suffer from 
unexpected and unexplained large finite size effects which need further study. 
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1. Introduction 

In the theory of first order phase transitions, a quantity of major interest is the 
interface tension between phases in metastable coexistence. In particular, in the study 
of nucleation phenomena one faces the problem to determine the curvature dependence 
of interface free energy/entropy between a nucleating droplet or bubble of the emerging 
stable and its surrounding metastable phase [H [21 [3l HI |5l [6l [7]. While an exact 
understanding of metastability in infinite systems is currently still lacking from the 
point of view of rigorous statistical mechanics [H |9], one can nevertheless try to obtain 
valuable information from the study of phase separation in finite systems [IDl [HI [12] • As 
analytical calculations are quite difficult for any nontrivial Hamiltonian, one resorts to 
simulations. For systems like Ising type spin models (lattice gases) [13] , binary mixtures 
[11] and simple fluids [151 [S] ? that have been studied by simulation, a convenient scalar 
order parameter like the total magnetization or particle number, which is extensive and 
even additive under partitions of the total system volume into subvolumes, is available, 
and whose density serves to distinguish between the different phases. Phase coexistence 
in an equilibrated finite system is characterized by the identity of the corresponding 
conjugate intensive quantity, whose physical meaning is that of an applied magnetic field 
or chemical potential. Thus one can extract the order parameter bulk densities of the 
coexisting phases from analyzing free energies obtained from Monte Carlo [IDl El [ISl [Tl] 
or molecular dynamics [161 [13 [H] simulations. In particular, armed with an additive 
order parameter, one can determine the equimolar volumes of the coexisting phases 
by evaluating the condition of vanishing adsorption of this order parameter for the 
corresponding choice of dividing surface. From such data it is straightforward to 
compute the interface tension [151 [H] • 

However, problems arise in cases for which a scalar additive order parameter is 
not available. On the one hand, the different phases may be characterized by different 
values of a multi-component rather than a simple scalar quantity (e.g. the formal vector 
of particle numbers of different species in multi-component fluids). On the other hand, 
there are cases where an explicit microscopic expression for an order parameter is frankly 
not known like e.g. in the study of protein folding. 

An extensive parameter that is always on hand and whose density - except for 
somewhat degenerate cases like hard spheres - may be utilized to distinguish different 
coexisting metastable and stable phases is the energy E. But unlike in the case of 
magnetization or particle number, the notion of an "equi-energetic" surface seems to 
be completely counter-intuitive, as attributing zero energy to the interface is in conflict 
with the fundamental physical principle that interfaces are regions with energy densities 
that are usually considerably larger than bulk ones. Thus, it is not clear how a 
reasonable measure of the corresponding subvolumes of the phases can be obtained 
from a microcanonical analysis based solely on monitoring the energy density. The 
purpose of the present paper is to show how this can be done, and that energy may still 
be a "good" parameter to determine interface tensions. 
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2. Review of the Grand Canonical Route to the Interface Tension 

We start by reviewing the traditional Gibbs dividing surface approach to the description 
of phase coexistence of two phases labelled a and /3 in a fluid |T9l [201 El 1211 [22] • Suppose 
flrst that the fluid has only a single chemical component. In spite of the fact that the 
problem of how to rigorously deflne coexistence of a droplet of phase a in (unstable) 
equilibrium with the surrounding /3 phase has not been solved up to date for the case 
of an inflnite system, it makes sense to consider such a situation in a large but finite 
system, in which such an equilibrium may actually be stabilized by imposing appropriate 
thermodynamic boundary conditions. Phase separation into two connected regions is 
then usually detected by observing regions of different density by monitoring the average 
inhomogeneous density proflle p(x). Following Gibbs one may then choose an arbitrary 
dividing surface, deflned as a levelling surface of zero thickness normal to the density 
gradient field, thus splitting the total volume V into subvolumes 

V = V^ + Vp. (1) 

Once this separation is agreed upon, any other extensive observable M can be split into 
homogeneous and so-called excess contributions as follows. If M assumes homogeneous 
equilibrium densities ma^mp in the phases a;,/9, we set 

M = + Mp + M^ (2) 

where 

Mo, = Vcjno., Mfi = Vfimp. (3) 

In a similar way, we construct the excess energy i?^, entropy S*^, Helmholtz free energy 
F^, grand potential fi^, particle number A^^', and so on from their homogeneous densities 
and the division ([T|). In a box of volume V = L"^, planar interfaces will usually form 
parallel to one pair of limiting walls. The excess quantities defined above will then 
generally depend on the chosen position of the dividing surface with area A = L'^~^. 
A notable exception is fi^, as can be understood from the fact that the corresponding 
homogeneous densities are the negative pressures of both phases, which in case of a 
planar interface must agree by simple stability arguments. Thus, for a planar phase 
separation geometry, the interface tension 

a = Q^A (4) 

is well defined, regardless of any parallel shift of the dividing surface. In contrast, the 
adsorption 

r = N^A (5) 

changes with such a parallel shift of A. In turn, the condition of vanishing adsorption 
r = then uniquely fixes the position of A and leads to an intuitively appealing 
definition of the "actual" position of the interface, known as the equimolar surface. 

While it is straightforward to generalize the definition of the position of such an 
equimolar dividing surface from the planar to that of a spherical or otherwise curved 
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case, the definition of tlie corresponding interface tension now requires considerably 
more care. For a spherical interface in = 3 let us agree to use a to label the phase 
inside the spherical volume and /3 to label the surrounding one outside the sphere. To 
stabilize a curved interface, the pressures on both sides of the interface must necessarily 
be different. In classical macroscopic physics this fact is encoded in the Laplace- Young 
(LY) equation 

Ap = po, - = ^, (6) 

which was derived from a mechanical analysis of the surface tension at the beginning 
of the 19th century. When promoting the statistical mechanics definition (j4]) from the 
planar case to that of curved interfaces, since Ap ^ for pi = Qi/V =: —Ui,i = 
the grand potentials densities cJq,,^/? of both phases will also disagree and so 

= n''{R) = Vuj -Va{R)iOa-Vf^{R)ujp (7) 

is bound to pick up a dependence on the radius. Thus, for a curved interface the interface 
tension a = o"(-R) appearing in ([6]) will itself be i?-dependent. At this point, however, 
it is crucial to realize that within the Gibbs dividing surface approach the choice of 
the radius R, being a purely theoretical construct, is in principle arbitrary, such that 
the classical relation ([6]) can only hold for a distinguished value of R. Nevertheless, 
physically observable quantities should not depend on the position of the artificially 
introduced dividing surface. Indeed, differentiating the definition of the excess grand 
potential in d dimensions with respect to R while keeping all other variables fixed, which 
is called a notional derivative and indicated by a bracket notation [d/dR], one arrives 
at the generalized LY equation [191 120] 

'da{R) 



Ap=(d-1)^ + 



dR 



(8) 



The classical LY equation ([6]) may only be recovered from this equation for the special 
choice of radius i? = i?,, for which 



da{R) 



^ 0. (9) 

R=Rs 



dR 

To account for this fact, the corresponding dividing surface is known as the surface of 
tension. Thus, Rs is a stationary point of C3"(i?) under a notional variation of R. In fact, 
it is not hard to see [T9l [20] that Rg is indeed a minimum of cr{R), which can explicitly 
be deduced from the general form of cr{R) in 2 or 3 dimensions 

according to which o"(-R) is universally determined from knowledge of Rg and cr^ = (j{Rs). 
The difference 6{Rs) := Re — Rs between radius of the equimolar surface Re and Rs, 
which has become famous under the name Tolman length, is known to be of molecular 
sizes. Of course, in principle one could also compute physical observables from e.g. the 
equimolar interface tension (Xe = cr{Re) or any other choice of radius R, since all physical 
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information is encoded in any such choice. However, the choice Rg is particularly 
convenient, as the condition (j9]) allows to condense many formulas to a considerable 
more compact and manageable form. To some extent this is also true for the choice 
R = Re, but the definition Tg = T{Re) = explicitly makes use of the particle number 
as an extensive parameter. For a multi-component fluid of r > 1 different chemical 
components, which we may label by an index t = 1, . . . , r, this complicates manners in 
a considerable way. In fact, the particle number N, adsorption F, associated chemical 
potential fi are then both promoted from scalar quantities to formal vectors N, F, fi, 
where N = (A^^, . . . , A^^) and so on. One now must deal with r-component averaged 
density profiles p(x) = (p^(x), . . . , p'^(x)). In general, such systems have complex phase 
diagrams, and if we concentrate on a particular transition, the levelling surfaces of 
different density components p*(x) may yield different equimolar surfaces for each chosen 
component, relative to which the remaining components form inhomogeneous adsorption 
layers. In principle it may still be possible to eliminate clumsy adsorption terms from 
formulas to some extent by an ingenious choice of dividing surface, the so-called Koenig 
dividing surface [23]. However, to fulfil its defining condition, all r components of 
the adsorption have to be balanced simultaneously with the multiple components of 
the associated chemical potential, which leads to profound numerical difficulties in the 
evaluation of simulation results. 

To illustrate these calamities, let us review the practical steps to calculate the 
interface tension from Monte Carlo free energy simulation data for a one component fiuid 
and compare them to the expected effort for a multi-component one. We choose a cubic 
simulation box of size N = L'^ with periodic boundary conditions. Carrying out grand 
canonical Monte Carlo simulations at some fixed temperature Tq chosen somewhat lower 
than the critical temperature T^, one first has to determine the coexistence chemical 
potential fi = yUo(^o)- Practically, this is done by implementing the equal-weight rule 
[21] numerically by performing a histogram re-weighting [25] to the simulated grand 
canonical probability distribution PtqVijl{^) of particle numbers and extrapolating the 
result to the thermodynamic limit L — oo. For a multiple component fiuid, these steps 
are already quite involved, as the approximately Gaussian-shaped peaks of Pj[)y^(N) 
are defined on a multidimensional space of variables and in presence of q possible low 
temperature phases one has to deal with g -|- 1 of them in order to pin down the r 
individual components of the vectorial chemical potential n to their coexistence values 



Nevertheless, suppose that this task has been carried out successfully. In the 
single component case, one next needs to resolve the detailed structure of PtoVmoI^) 
by e.g. Wang-Landau sampling, eliminating residual inaccuracies by a weighted Monte 
Carlo production run. In this way, one obtains a dimensionless excess free energy density 



which depends on the scalar density p = N/V (we have dropped a purely To-dependent 
normalization part which is irrelevant for what follows; the usual thermodynamic 
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notation Pq = I/Aj^Tq should not lead to any confusion with the label P for one of 
the two phases). At a first order phase transition, the finite size potential (3Qf^^\TQ,p) 
has a double-welled shape with a flat central plateau, and an analysis of its fine details 
reveals several distinct ranges of the total density p between its two minima, for which 
one may observe phase separation in planar, cylindrical (in d = 3) and spherical shapes 
induced by the imposed periodic boundary conditions. These regions are also detectable 
in the derived finite size canonical "excess" chemical potential 



ii"''(n,p)^ . (12) 



V 9p 



To 



From the distorted "Van der Waals loop" shape of p^^^Tq, p), one can infer that for 
small enough p the equation p^^\Tq, p) = p generally allows for at least three roots 
p^a\TQ, p) < p^^\Tq, p) < p^^\Tq, p) corresponding to two bulk densities Pa\ p^p'^ of 
coexisting phases a, /3 at total density p. This allows to introduce the total dimensionless 
grand potential density 

a;(^)(To,/i) = /^^)(To,p(To,/i)) -/ip(^)(To,/i) (13) 

and the grand potential densities 

ujf\T,,p) = f''^\T,,p^{T,,p))-pp\t\n.p), (14) 
ujf\T,,p) = P\To,p^iTo,p))~ppf\To,fi), (15) 

from which it is straightforward to compute for a given volume partition V = 
V^(i?)+Vfl(-R), and thus, using (jlj), the notionally i?-dependent interface tension a^^\R). 
In detail, Rg and ai^^ are found by numerical minimization of o'{R) with respect to R, 
while ai^^ can be calculated from the spherical equimolar volume determined by the 
lever rule 

VajRe) ^ P- Pa V^jRe) ^ Pfi - P ^^^-^ 

V Pl3- Pa V Pli- Pa 

for Va{R) = 4:71 R^/3 in d = 3 and Va{R) = nR^ in d = 2, respectively. In passing we 
note that since ^l^{Re) = F^{Re), ai^^ is exclusively determined by /*^^^(To, p(To, /i)) 
(cf. [151 E]). 

All this is very fine - however, in a multi-component setting the numerical effort to 
carry out these steps successfully is again forbidding. What is needed is an intrinsically 
scalar approach. 



3. Microcanonical Approach 

As mentioned above, a basic extensive scalar quantity that is always available is the total 
energy E of the system. It is by now well known that while phase transitions are only 
well defined in an infinite system in the strict sense, they can nevertheless be studied 
conveniently by analyzing the so-called "convex intruder" in the microcanonical entropy 
S{E) of a finite version of the system [261 EZj- By its very definition, such a convex 
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E, 



E 



Figure 1. Schematic illustration of a composite entropy with a convex intruder. 



intruder in S{E) gives rise to a corresponding anomaly of the microcanonical inverse 
temperature (3{E) and a possibly negative branch of the accompanying microcanonical 
specific heat with quasi-singularities at the intruder boundaries (for a schematic picture, 
see Figure [1] below) . In a finite (or at least in some sense "small" [26] ) system such an 
observation does not contradict thermodynamic stability requirements [28]. One now 
observes that up to a sign, the overall appearance of (3{E) is quite similar to that of 
the excess chemical potential jl^^\(3^p) discussed above. Indeed, parallels between the 
former constructions and the ones carried below are not accidental. 

Recording the energy probability distribution PtqVij.{E) at fixed Tq for different 
values of /i, we can locate the coexistence chemical potential /^q = A*o(^o) by applying 
the equal weight rule [291 1211 EQ] to the two separate approximately Gaussian-shaped 
peaks, as will be discussed below for the case of the Potts model. In any case, this 
may be feasible even in case of a multi-component order parameter, since the domain 
of PtqVii{E) is the scalar quantity E. For the rest of the simulation, the parameters 
To and /^g remain fixed. It is then convenient to formally replace the energy E hy a, 
grand-canonical version 



in the same way as one may formally replace the original Hamiltonian by a grand- 
canonical one in computing the grand-canonical partition function. For our purposes, 
the "background" homogeneous chemical potential /Zq can thus be eliminated from the 
list of thermodynamic variables. With /Xq fixed, we now record the number of states 



with grand canonical energy E m. a, fiat histogram simulation followed by a weighted 
Monte Carlo production run. Once g{E, V) is known, we arrive at the grand canonical 
partition function in the form 



E-^IqN -> E 



(17) 



g{E,V)^6 



,s(Ey) 





(19) 



E 



e 
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where we introduced the energy and entropy densities e = E/V and s^^\e) = 
S{E,V)/V , vahd for arbitrary /3. Suppose that at a particular inverse temperatures 
(3 the sum is dominated by its largest summand in the limit of large system size. The 
corresponding energy density e(/3) is determined by the equation 



13 



V de 



V 



/3(^)(e)L=,(.)(,^, (20) 



=e(i')(/3)) 
E(i){/3) 



which expresses the equality of inverse canonical temperature /3 and microcanonical 
temperature (3^^\e) for this special value of e. One may then approximate 
PQ{P,V,fi,) = - \ogZ{P,V,^io) 

/3a;W(/3, /Xq) ^ /^e^^H/?) - ^^^He^^H/?)), (21) 

i.e. as a Legendre transform. This equation is quite similar to f lT3|) . — ^'•^•'(e) playing 
the role of f^^\ e that of p and —/3 that of /i. Similar to the strategy employed above 
for the chemical potential, for a certain range of inverse temperatures /3 the equation 
/3 = /3^^\e) has (at least) three different solutions Ca < e < corresponding to the total 
energy density e and those of the two coexisting phases a, (3. If we now let (3 approach 
this interval around /3q, the presence of the convex intruder in s^^^(e) makes the above 
discrete saddle point approximation break down since not one but at least three "saddle 
points" are found, which in the non-degenerate case of three roots correspond to the 
dimensionless grand potential densities 

/3a;(^)(/3)=/3e-.(^)(e), (22) 
/3a;(^)(/3)=/3e,-s(^)(e„), (23) 
f3uf\p)=Pep-s'^'^\ep). (24) 

At this stage, we have reached our goal announced above, as these quantities, once 
they are determined, allow to compute the excess grand potential and thus the interface 
tension a^^\R) from formula (171). By minimization of a^^\R) w.r.t R we obtain the 
radius of the surface of tension Ri^^ and the corresponding surface tension ai^\ A 
calculation of Re is, of course, not feasible in this approach. 



4. The 2d Potts Model 



We illustrate our strategy taking the example of the g = 30 nearest neighbour Potts 
model in c? = 2. The g-state Potts model [31], whose Hamiltonian on a simple cubic 2d 
lattice oi N = L"^ sites in zero external field is given by 

?/[{.(x)}] = 5^[l-5,(.),,(y)], (25) 
(xy> 

where s(x) G {1, . . . , g} is ideally suited for this purpose for a number of reasons, (i) 
First of all, for g > 4 the model undergoes a temperature-driven first order phase 
transition, (ii) Regarding this transition, a wealth of rigorous results [311 1321 l33l [311 (35] 
is available in the literature which can serve to benchmark our simulation results. For 
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the first order phase transition temperature of a bulk system, one has the exact analytic 
expression 

l/To = /3o = ln(l + v/g). (26) 

In [361 El] it was reported that, as was expected from general arguments, the inverse 
temperature 1/Tq{L) = Pq{L) at which the ratio of the two weights of the thermal energy 
probability distribution is just g, agrees with the exact bulk value (5q = (3o{L = oo) up to 
exponentially small corrections. Thus, Po{L) serves as a convenient definition of a finite 
size transition temperature. Other rigorous results include the latent heat per volume 
[33] , the limiting internal energy densities at Tq and the difference of specific heats |32j . 
Furthermore, the reduced interface tension (i.e. the interface free energy density at (3q^ 
multiplied by (3o) between the disordered and one of the ordered phases along the square 
lattice (10) direction was rigorously determined [35] to be 



where 



with 



2o-o/d 



n=0 



In 



(29) 



Equally important for us is the fact that there is even a rigorous calculation of the full 
anisotropic interface tension available [37]. However, since these calculations are too 
involved to be reproduced here, we content ourselves with noting that at /3q the resulting 
anisotropy for g = 30 calculated from the formulae in [37] is vanishingly small. This 
happenstance is a very important prerequisite for any attempt to apply our evaluation 
strategy for the interface tension, which rests on a presupposed spherical symmetry of 
bubbles and droplets. 

(iii) The q state Potts model's order parameter is not scalar, but has dimension 
g — 1. Since this is a central issue in the present context, let us briefly review its nuts 
and bolts. Guided by physical intuition, a scalar "order parameter" m could be defined 
by the following reasoning [an |29]. Let N^°-^ denote the number of spins of a given 
microstate with value s(x) = a, where 1 < a < g. Let A^'max := ma,x{N^^\ . . . , N^'^^). 
Then 

M := (30) 
g- 1 

Obviously m := M/N is confined to values < m < 1, and m = for complete 
disorder, while m = 1 for any perfectly ordered domain. In the Ising case q = 2, m 
indeed corresponds to the modulus of the magnetization density of the system. Thus, 
if we break up the system volume into subvolumes Vi and add up their different Mj's, 
generally M ^ ^ . Mj, i.e. M is not additive between subsystems. 
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On the other hand, in Zia and Wallace construct a full (g — l)-component order 
parameter. They introduce q unit vectors in e^"^ G a = 1, . . . , q, q > 1, such that 

the following relations are satisfied: 

eWeW = ^^3^. (31) 

Any set of such vectors defines a generalized tetrahedron in M''"^, the q vectors pointing 
from the center to each corner. Trivially, the vectors e^"^ cannot be linearly independent. 
Instead, they satisfy the geometrically evident sum rule 

^e(")=0. (32) 

a 

With the one-to-one correspondence s(x) e^^^^^'^ =: e(x) understood, one associates 
a g — 1 component order parameter 

M:=^e(x) eM."-^ (33) 
xer 

with each given microstate, which we will call the magnetization. In terms of the 
occupation numbers A^'-'*^ of the Potts spin states 

Q 

M = ^A^^'^^e^'^). (34) 

a=l 

If one multiplies (!34|) with any of the unit vectors e^''^ it is easy to see that 

M") := e(')M = (35) 

q-1 

Clearly, m(^) := M^/A^ G [-l/(g-l),l] agrees with m as defined in (l30l) provided 
= A'^tnax- This clarifies the role of m as well as the low temperature domain structure 
of the model. Namely, suppose that AtC") = N^^^. Then, we can rewrite f l34D as 

M = V(A^(^) - Ar("))(-e('^)). (36) 
' 

The q domains D^''^ are thus geometrically represented by convex cones enclosed by the 
set of vectors {(— e'^'^^) : a = 1, . . . , q, a ^ b}, and (l35l) gives the projection of M onto 
the average direction Yla^b('~^^'^^) ~ ^^'^^^ symmetric group of permutations of q 
numbers acts as the underlying symmetry. 

At small values of m, several very small fluctuation "clusters" of the same or 
competing spin values may coexist, causing the system to jump randomly from one 
domain cone to another. On the other hand, the parameter m^^^ is always additive by 
construction, but does depend on the chosen direction b in M-space. Only for values 
of m^''^ larger than a certain threshold do we find agreement of the order parameters 
computed from ( l30l) and ( l35ll . since then the direction along which the projection from 
M to M occurs is uniquely determined by the value b of the majority occupation 
number. A cluster decomposition performed during the course of a simulation can 
provide information in analyzing these fluctuations in order parameter topology. 
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Thus, the parameter (15U]) is only additive for magnetizations Mi and M2 both 
belong to a single common domain D^'^\ and thus cannot be used meaningfully as a 
parameter in a Gibbs dividing surface construction. On the other hand, sampling the 
free energy as a function of the full order parameter M is out of the question. In 
other words, we are exactly in the situation outlined earlier, and thus embark on a 
microcanonical strategy instead. Well, not quite. Actually, the situation for the Potts 
model is not as complicated as the one outlined in the introductory section. This 
is largely due to the fact that there is no need to determine a coexistence "chemical 
potential", which would correspond to an external vectorial magnetic field. As the 
different g-spins of the Potts model are all coupled in the same way, this external field 
is fixed to be exactly zero. 




Figure 2. q = 30 Potts model at in 2 dimensions with periodic bomidary conditions at 
L = 100: snapshot of typical "droplet" configuration at scalar order parameter value 
m « 0.06. 

5. The Microcanonical Heat Bath Algorithm 

Any successful design for a Monte Carlo algorithm devoted to the study of phase 
separation at phase transitions of pronounced first order character must address the 
phenomenon of exponentially diverging autocorrelation times 

r(L,/3) ~exp(2/3crooi:'-') (37) 

accompanying first order phase transitions which known as hyper- critical slowing down 
[3H] (in fl37p a d- dimensional cubic box of volume L'^ with periodic boundary conditions 
is assumed). A Wang-Landau type of algorithm is capable of overcoming large entropy 
or free energy barriers separating different stable or metastable phases [ID] and is in 
principle straightforward to implement. In its original microcanonical version, the 
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algorithm directly yields the density of states (or, for a discrete system, rather number 
of states) 

g{E) = J2S{E.-E) (38) 

V 

of a system having microstates v with energies (we are somewhat casual about 
the use of the Dirac delta function for simplicity, and we have put /cs = 1 for the 
same reason) with high precision, which is related to the microcanonical entropy by 
SiE) = In g{E). To control possible residual errors, one may thus determine an 
approximate microcanonical density of states by Wang-Landau simulations and perform 
subsequent biased Monte Carlo production runs with statistical weights based on 
this approximate density of states. In fact, knowledge of S{E) conveniently allows 
to determine a multitude of other T-dependent observables at various temperatures 
simultaneously with high precision, as is explained in more detail below. 

It remains to construct a suitable move set for a microcanonical Wang-Landau 
simulation scheme. Single g-spin updates are simple to implement and may be a 
reasonable choice for Ising systems far from criticality, but are quite inefficient in 
exploring the regions of phase space of the large q Potts model which are of interest 
for studying phase separation, namely those configurations, in which a single ordered 
domain of Potts spins s(x) of, say, value s(x) = q coexists with a disordered background. 
Indeed, suppose that during the course of the simulation, our random walk arrives at a 
particular configuration, in which almost all Potts spins agree with this condition, while 
just a few, say, s(xj), i = 1, . . . , k are yet disordered. In such a microstate, chances are 
only A^~^ X {q — 2)~^ that of the "missing" sites s(xj) is indeed drawn and its spin s(xj) 
be assigned the "right" value q in creating the next trial configuration, thus matching 
the surrounding domain. 
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Within the canonical ensemble, it is well known that the heat bath algorithm ^T[W2\ 
is superior to the standard Metropolis scheme for high q Potts models. In detail, let 
^(m) denote the total energy of the system in the microstate fx. Choose a random site 
X and let := s(x) denote the value of the Potts spin variable at this particular site. 
Furthermore, let yi,i = 1, . . . ,z denote the nearest neighbour sites of x. Defining the 
local energy at site x by 

z 

Eq = -E'(iocai)(g, := - X] ^q,s{yi}^ (^9) 

one can split 

i?('^^=i?(local)(gM,x)+i?rSt (40) 

and define a set of heat bath probabilities 

which are manifestly independent of the value of the initial central Potts spin s(x). The 
canonical heat bath algorithm amounts to choosing a new value q^, G {0, . . . ,g — 1} 
for this spin with probability in every step. It is easy to see that the resulting 
algorithm satisfies detailed balance as well as ergodicity. In terms of generation and 
acceptance probabilities, we have g{fi = ^ and a(/i u) = 1^ i.e. the stochastic 
character only enters in the generation of configurations, which are, once generated, 
always accepted. 

It is straightforward to translate these ideas from the canonical to the 
microcanonical setting. To illustrate the correspondence, let us denote the canonical 
Boltzmann weights by 

TT, := ^— — . (42) 
Then the canonical heat bath probabilities Pq can quite trivially be rewritten as 

g— /3(_Bg+_ErCSt) yj-^ 

e-/3(S.+-E;rest) " V^-l • (43) 

Now, to translate the algorithm to the microcanonical ensemble, we simply replace 
VTj, — )■ 1/ g{Ey). In terms of the microcanonical entropy S{E) := In (?(£'), we can rewrite 
the above probabilities as 

(44) 

Once we have determined the microcanonical entropy S{E), it is in principle 
straightforward to obtain the free energy as a function of T and some other (preferable 
scalar) observable o = 0[{s{x.)}], where (9[{s(x)}] denotes e.g. the magnetization m, the 
projection of the order parameter m^""^ along an arbitrary fixed internal direction a, the 
size of the largest geometric or Swendsen-Wang cluster, all simultaneously computed for 
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different temperatures from the underlying microstates {s(x)} visited during the course 
of a single microcanonical biased Monte Carlo simulation. 

To sketch this procedure, we consider the constrained microcanonical density of 
states g{E, a), which is formally written as 

g{E, a) =J2Ho- 0[{s{^)}]) 6 {E - ?/[{.(x)}]) . (45) 

{^(x)} 

The corresponding (conditional) probability to find the value o of (!?[{s(x)}] at total 
energy E is 

where h{o, E) denotes a two-dimensional histogram recorded during the course of the 
simulation. But, according to the rules of conditional probabilities, this precisely implies 
that the canonical probability to find the value o of C[{s(x)}] at inverse temperature 

13 is 

P(o|/3) = P{o\E)P{E\P) oc h{o, E)P{E\P). (47) 

E E 

To obtain this probability at any given temperature from a microcanonical Monte Carlo 
simulation biased by the predetermined density of states g{E), we thus only need to 
reweight the recorded histograms h{o,E) by the known function P{E\I3). The desired 
constrained free energy density 

f{(3,o) = -J^lnJ]5(o-0[{.(x)}])e-^«[i^W>l (48) 

^ {^(x)} 

can be recast in a similar way, since 

^5(o-0[{s(x)}])e"^«K^W>l 

{^(x)} 

= E^"'"" J2^io-0[{si^)}])6iE-n[{si^m 

E {.(x)} 

= Z{P)Y,P{E\P)g{E, a) = Z(/3)P(o|/3), (49) 

E 

and is thus (up to an unimportant constant) given by 

/(Ao) = -^lnP(o|/3). (50) 



6. Microcanonical Results 



We have conducted a series of Landau- Wang simulations followed by weighted Monte 
Carlo production runs of a 2d square lattice Potts model of size N = L"^ to determine 
the density of states g^^\E) and thus the entropy density s'^^\e) = N~^\ng^^\E/N) 
and the microcanonical temperature f3^^\e) = ds^^\e)/de. Since we are interested 
in the development of phase separation in small to finite system sizes, we carried 
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Figure 4. Microcanonical entropy densities s^^^(e) for the q = 30 Potts model in 
d = 2 dimensions with periodic boundary conditions for various (but indistinguishable) 
system sizes. 
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Figure 5. /3^^'(e) with energy density ranges of droplets and bubbles marked; the 
horizontal line displays /3o ~ 0.86829. 

out simulations for linear sizes L = 40, 50, . . . 100. For the g = 10 Potts model, the 
resulting signs of phase separation were not observed to be very pronounced. However, 
increasing g to g = 30 clearly revealed the expected convex intruder. But even then, 
from merely looking at the entropy density s^^\e) it is virtually impossible to detect the 
delicate features appearing at finite system sizes that we are interested in (cf. Figure H]). 
However, numerically taking the derivative of s^^\e) with respect to e, we obtain the 
microcanonical inverse temperature /3^^^(e) which provides a detailed view of the delicate 
substructures hidden in s^^\e) (cf. Figure [5]). Inspection of the resulting curves gives 
a first hint on the quality of our simulation data, as one should take into account that 
numerically differentiating potentially noisy data should greatly magnify any statistical 
irregularities and errors in such data. To resolve the convex intruder in the original 
entropy data s^^\e), it turns out to be convenient to consider the auxiliary dimensionless 
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function [13] 

A(^)(E,/3)/iV = A(^He,/3) := /3e - s(^)(e). (51) 

Of course, A'^^^(e,/3) coincides up to a constant with the logarithm of the canonical 
energy probability function P^^\E\/3) at inverse temperature (3. However, for our 
present purposes we may regard A*-^^ as a finite size "Landau potential", i.e. an 
incomplete Legendre transform of the microcanonical entropy, and compare its features 
to those of canonical Landau potentials (cf. [HlllS]). Thus, let us tune the parameter 
P to values near the bulk inverse transition temperature Po and analyze the resulting 
shape of X^^\e, (3) as a function of e. As expected, for such temperatures \^^\e, f3) 
resembles a somewhat distorted double-well shape (cf . Figure [6]) with two pronounced 
minima at energy densities Cc = E^/N and = E^/N (the subscripts "c" and "v" 
correspond to "condensed" and "vapor") separated by a large "Landau free energy 
barrier" with a "flat" central plateau of practically constant and vanishing or at least 
quite small slope. We identify Cc and with the equilibrium energy densities of the bulk 
"condensed" (ordered) and "vapor" (disordered) phases, while the thermodynamics of 
their possible coexistence configurations in encoded in the features of the potential well 
between them. The fiat central region of the potential, which signals phase separated 
configurations with a slab-like interface geometry [H], is, of course, also reflected in the 
central linear section of I3'^^\e) at the level of I3'^^\e) ~ (3q (cf. Figure [5]). 

Apparently, Figure [6] also illustrates the dilemma of defining a "proper" finite 
size transition temperature for a system with a highly degenerate low energy domain 
structure. On the one hand, one could naively try to adjust /3 to such a value that 
both minima of A*-^^(e,/3) are of equal height. This choice precisely corresponds to 
the "equal height rule" for P^^\E\T). However, at such a temperature, one observes 
a noticeable slope in the central "fiat region" of \^^\e^f3)^ which signals that phase 
coexistence is not well established. In a plot of the quasi-Gaussian function P^^\E\T) 
this and other delicate features outside the peak regions do not give themselves away 
to the naked eye, since they are exponentially suppressed. On the other hand, choosing 
the inverse temperature /3 to agree with the rat io-of- weights temperature with L — )■ oo, 
which, as discussed above, converges exponentially fast to the exactly known inverse 
bulk transition temperature [33] 

/3o = ln(l + ^) ^'=°' 1.86829, (52) 

the fiat central region of A''^''(e, /3o) is found within numerical precision to be horizontal, 
i.e. with vanishing slope, but now one notices a pronounced difference in height between 
the two minima at energy densities ec,e„, which diminishes with growing system size. 
Numerically, this height difference is seen to approach 

A(^ne.,/3o) - A(^nec,/3o) ~ iV"Mng. (53) 

Physically we can interpret these findings as follows. Suppose that precisely at the 
inverse transition temperature /3o the system initially starts out in an ordered equilibrium 
state. Then, the probability to generate a fluctuation yielding a mixed state where half 
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Figure 6. A'-^^(e, /3o) with energy density ranges of droplets and bubbles marked 



of the available volume is turned into a disordered state separated from the ordered part 
by a straight interfacial line, should differ from the "inverse" probability to produce the 
same state starting from the disordered state by a factor of q, since in the latter case 
q possible ordered states are available, while only one disordered configuration can be 
formed in the former one. Taking the logarithm of the fraction of these probabilities 
then produces (!53|) . This completes the picture, since, as could have been anticipated 
from the practically perfect Gaussian nature of the two peaks in P{E\1/ Pq), the ratio- 
of-weights and ratio-of-heights temperatures are found to numerically agree for practical 
purposes. Without going into the details we also note that the equal height temperature 
can also be shown to coincide with the inverse temperature found by imposing a 
microcanonical Maxwell construction, i.e. chosing the value of /3 for which suitable 
defined areas obtained from integrating P{E) between Ec and E^ coincide |16]. 

At this point, a few additional comments concerning the true nature of the 
above "Landau potential", the non-monotonous behaviour of (3^^\e) and the resulting 
appearance of branches of "negative specific heat" are in order. In fact, in the 
thermodynamic limit, /3(e) = limi_i.oo /3^^'*(e) indeed decreases with e up to e = Cc, 
stays constant at /3(e) = /3o up to e = e^,, and then decreases further, as it should be: no 
trace of any metastable states (or even "unstable states") is left in the /3(e) curve. Of 
course, this must be so: apart from statistical errors, Monte Carlo simulations yield the 
equilibrium statistical mechanics of any such model Hamiltonian exactly, and metastable 
or unstable curves cannot be the output of exact calculations in the framework of 
equilibrium statistical mechanics. So for L — )■ oo the minimum position e^jl^ of /3*^^)(e) 
moves towards Cc (and its depth vanishes); similarly, the maximum position emix moves 
towards e^ (and its height, relative to /3o, vanishes as well). It is interesting to recall 
the physical significance of these extrema: for ec < e < e[^^^ the finite system is still 
homogeneous, and the minimum is the signature of the first appearance of a "bubble" 
of the disordered phase within the otherwise homogeneously ordered phase (cf. Figure 
[3]), while the maximum is the signature of the first appearance of a "droplet" of the 
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ordered phase within the otherwise homogeneous disordered phase (Figure [2]). As long 
as such "heterophase fluctuations" are absent, finite size effects are small in the curve 
P^^^ (e) shown in Figure[5l the strong finite size effects in between e^^j^ and Cmlx are due to 
interfacial contributions to the "Landau potential" (Figure[6]), which are of relative order 
1/L in Figure El Thus, the branch of negative "specific heat" resulting from Figure [5] in 
the region where (3^^\e) is an increasing function of e is not at all an unphysical result, 
but simply refiects the importance of interfacial free energies in finite microcanonical 
systems. In addition, our use of the nomenclature "Landau potential" merely refers to 
the incomplete character of the Legendre transform (ISTI) . but should not mislead the 
reader to confuse this potential with "Landau potentials" of similar appearance as they 
are constructed in mean-field theory. In fact, our potential A^^-*(e, (3), whose information 
content is, after all, identical to that of the full microcanonical entropy density s^^\e), 
describes the thermodynamics of two-phase coexistence in an inhomogeneous finite 
system without any approximation, and thus is conceptually quite different from a mean- 
field potential, which is constructed under the implicit constraint that the system is in 
a homogeneous phase throughout, whereas such states are thermodynamically unstable 
in reality. 

We can gain confidence in the overall correctness and general quality of our data 
by comparing the exact value = 0.29277 of the reduced g = 30 planar interface 
tension as calculated from (1271) to the one obtained by a finite size extrapolation of 
our data. In fact, as there are two minima of A*^^^(e,/3o) at energy densities ei^\ei^^ 
whose values differ by ~ In30/A^ as discussed above, there are two corresponding sets 
of data {A(^)(eL''ix, /3o) - A(^)(ei''\ /^o)}, {A(^)(eL'^i., /3o) - A(^)(e^''\ /3o)}, corresponding to 
the difference between the central barrier {\^^^ [cmlx, Po)} taken at some energy density 
e^mL ^ {ei^^ + ei^^)/2 and the left and right minima {A(^)(eF, /3o)}, {A(^)(ei^\ /3o)}, 
respectively (see Figure [6]). A standard finite size scaling extrapolation of these data to 
L — 7- oo in the form 

(3o(^c,v{^) = /3oO-c,i,(oo) + const/L, (54) 

which is displayed in Figure [71 gives the two values Poadoo) = 0.292168 and 
/3oO"t,(oo) = 0.291441 for the left and right difference, respectively, whose average 
/3o(crc(oo) + crt,(oo))/2 = 0.291805 differs by less than 0.4% from the exact value 
PoO'oo = 0.29276 computed from Formula ( 127]) . 

From (jl]), ([9]) and (!22|) - (!2^ it is now straightforward to compute a^^\Rs) and 
Rs for any prescribed total energy density e. Note, however, that these formulae are 
only valid for a spherical phase separation geometry. The corresponding approximate 
density regions within which one may expect states which on average resemble spherical 
droplets and bubbles to dominate in the sampled microscopic system configurations may 
be found by visually inspecting the slopes of A*^^^(e,/3o) and (3^^\e) in Figures [5] and 
El respectively, and cross-checking these ranges by examining corresponding snapshots 
taken during the course of the simulation. The total energy density regions for which 
we may expect the appearance of spherical droplets and bubbles are indicated in colour 
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Figure 7. Finite size extrapolation ([54|) of reduced interface tension for planar 
interface. 
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Figure 8. Normalized interface tension (t{Rs)/(Joo of droplets. The ranges of radii, 
for which the curves actually describe spherical droplets are marked in colour. 

in Figures |5] and El 

Our results for the interface tension a^^\Rs) at the surface of tension are gathered 
in Figures |8] and HI In correctly interpreting these results, it is quite important 
to understand that they have been obtained from (jlj), and ( !22l) -( IMl) under the 
assumption of spherical geometry. The corresponding ranges of inverse radii for which 
one can expect this assumption to be valid have been marked in colour in these figures. 
Outside of these ranges, the data do not accurately describe a physical interface tension, 
but merely serve as a guide to the eye. 

Looking at these results, one instantly notices the large finite size effects, 
manifesting themselves in the considerable offsets between the consecutive considered 
L- values, which strikingly fail to collapse onto a common "master curve". Currently, 
we find it difficult to understand the origin of this behaviour. For the planar interface 
tension, which is described by our data quite accurately as discussed above, strong but 
regular finite size effects are indeed expected. They may be heuristically understood 
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Figure 9. Normalized interface tension (j{Rs)/<Joo of bubbles. The ranges of radii, 
for wliieh the curves actually describe spherical bubbles are marked in colour 

in terms of the L-dependent truncation of the wave vector spectrum of capillary waves 
running parallel to the interface. However, for a spherical bubble or droplet, identical 
radii R should yield identical values of the interface tension, once the surrounding box 
has been chosen large enough to kill finite correlation length effects, which are however 
expected to be vanishingly small for a strong first order phase transition. 

At the moment we do not have a clear explanation of these strong finite size effects 
which prevent us from a further meaningful analysis of the curvature dependence of 
a{R). In similar approaches to study the interface tension of curved interfaces in a 
truncated 3d Lennard- Jones fluid [U] and a 3d fee lattice gas model |18], we also 
find certain finite size effects, but they are much less pronounced than those of the 
present case. However, we have also observed finite size effects of comparable size in 
computing the interface tension from canonical simulations of a 2d Ising model. Thus, 
we believe that the large magnitude of the finite size effects has nothing to do with 
our microcanonical approach, but is rather related to the fact that both systems are 
two-dimensional. With growing linear system size L, the gaps between consecutive 
a^^\Rs)-cnTves recorded in Figures |9] and |9] obviously diminish, so these curves are 
expected to eventually collapse onto a single "master curve" for large L and Rg. On the 
other hand, for L > 110 we report that our Monte Carlo production runs failed to be 
sufficiently ergodic, indicating large residual entropic barriers with respect to "hidden" 
observables beyond our one-dimensional energy-based sampling, while at the same time 
the barriers observed in T^^\E) had already risen to a value of 60. To study much larger 
systems would thus require to overcome these additional hidden barriers, presumably 
by employing much more elaborate sampling techniques than the ones we are using here 
(see e.g. |19] or [50] for promising approaches). However, in our current work we are 
not interested in exceedingly large droplet sizes and their accompanying huge entropy 
and free energy barriers. Rather, out intention is to focus on the behavior of droplet 
and bubbles of moderate size, as this is the only regime that is of practical relevance 
for nucleation related questions. In any case, the origin and nature of the encountered 
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finite size effects must currently be left to further study. 
7. Conclusions 

In this paper, we have addressed the investigation of phase coexistence of systems with 
a more-component parameter in the context of computer simulations, which necessarily 
involve systems of finite size. Such simulations of phase coexistence often are done 
with the motivation to extract information on the interfacial tension of fiat and curved 
interfaces. While for systems with a scalar (i.e. one-component) order parameter 
this problem is normally considered in the grand-canonical and canonical ensemble of 
statistical mechanics, we have given a concise discussion of this approach to show that its 
extension to the multi-component case is formally possible but practically unfeasible. 
We then have presented, as an alternative, a microcanonical approach based on the 
number of states g{E, V) of energy E for a system having a finite volume V. In the 
entropy versus energy curve S{E) for the finite system there is a convex intruder (Figure 
[I]), and the idea we follow in the present paper is to carefully analyze this intruder 
as a function of system size, in order to extract information on interfacial tensions. 
We exemplify our approach for the two-dimensional g-state Potts model with a large 
number of states {q = 30), proposing also an extension of the heat bath algorithm from 
the canonical to the microcanonical ensemble. From these simulations we obtain very 
precise information on S{E) and also an effective potential X^^^e, f3o) per lattice site, e 
being the energy density and /3o the inverse temperature where in the thermodynamic 
limit (V = — > oo) the first-order transition from the ordered phase to the disordered 
phase occurs (Figure |6]). Also the derivative (3^^\e) = d\^^\e)/de is obtained with 
meaningful accuracy (Figure |5]). We have shown that the loop in such l3^^\e) vs. e 
curves has nothing to do with the "van der Waals-like" loop of mean-field theories: 
in the latter, such loops describe a path of homogeneous states connecting the two 
phases between which the transition occurs; in reality, our loops (Figure [5]) reflect two 
phase coexistence in finite systems, all parts of the loop describe full stable thermal 
equilibrium; any interpretation in terms of metastable or unstable states would be 
completely misleading. There is nothing mysterious about the "negative specific heat" 
that often is attributed to such loops - the whole loop just reflects interfacial effects, 
just as the "hump" in between the two minima of the "Landau potential" in Figure El 
all these features disappear proportional to 1/L in the limit L — t- oo, and the correct 
horizontal parts in between Cc and e„ remain, as it should be. Thus, one should not 
be mislead by mean-field concepts when discussing first-order phase transitions in finite 
system in the microcanonical ensemble. 

We have found that in the fiat region in the center of Figure [6] the data allow 
an accurate estimation of the interfacial tension of fiat interfaces between ordered and 
disordered phase (Figure Ej), although also in this case finite size effects are clearly rather 
pronounced, and an extrapolation to L — )■ oo is mandatory. However, the analysis of 
the ascending parts of X^^\e, (3o) in Figure [6] in terms of the radius-dependent interface 
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tension of droplets (Figure |H]) and bubbles (Figure ^ is more subtle: again huge finite 
size effects occur, and it is not possible at fixed radius Rg to extrapolate to L — )■ oo, 
because due to the droplet (bubble) evaporation/condensation transitions droplets at 
fixed radius Rs are only stable in a rather restricted range of L. While naively one 
could expect that different choices of L yield mutually compatible results for a{Rs), 
as approximately happens for one-component systems in d = 3 dimensions, this is not 
the case here. Of course, our analysis does not explicitly consider the fact that at a 
given value of e and the corresponding average value of /3 (Figure |5]) in the two-phase 
coexistence region at a given value of L the droplet (or bubble) is strongly fluctuating 
both with respect to its size and its shape (Figures [2] and [3]). We assume the shape 
of the droplet or bubble to be spherical, otherwise the information recorded does not 
suffice to extract a{Rs). Future work along such lines must analyze this problem of 
droplet (bubble) fluctuations more closely, possibly by recording additional observables 
related to the droplet or bubble, to allow estimation of o"(i?s) within reasonable error 
limits. In fact, in d = 3 the fluctuations are found to be indeed much less pronounced, 
and - at last in the one-component case - meaningful results for cr{R)s are accessible 

[miiH]. 
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